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Abstract 



The statistical evidence (or marginal likelihood) is a key quantity in Bayesian statistics, allowing one to 
assess the probability of the data given the model under investigation. This paper focuses on refining the 
power posterior approach to improve estimation of the evidence. The power posterior method involves 
transitioning from the prior to the posterior by powering the likelihood by a temperature variable. In 

i ' 

d ' common with other tempering algorithms, the power posterior involves some degree of tuning, and this 

^ , ' paper addresses this issue. The main contributions of this article are twofold - we present a result from 

the numerical analysis literature which can reduce the bias in the estimate of the evidence by addressing 
the error arising from numerically integrating across the temperature. We also address the choice of 
ON • temperature ladder, and present an adaptive algorithm which gives excellent performance in the examples 

£f~) ' considered here. A key practical point is that both of these innovations incur virtually no extra cost. 

ON" Keywords: Marginal likelihood, Markov chain Monte Carlo, Power posteriors, Statistical evidence, 

£Nj . Tempering, Thermodynamic integration. 



^ ■ 1 Introduction 

The statistical evidence (sometimes called the marginal likelihood or integrated likelihood) is a vital quantity 
in Bayesian statistics for the comparison of models, mi, . . . , m\. Under the Bayesian paradigm we consider 
the posterior distribution 

p(9i,mi\y) <xp(y\9i,mi)p(6i\mi)p(mi), for i = 1, . . . , /, (1) 

for data y and parameters 9{ within model mi, where p(9i\m,j) denotes the prior distribution for parameters 
within model rrii and where p{rrii) denotes the prior model probability. The evidence for data y given model 
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mi arises as the normalising constant of the posterior distribution within model mi, 

p(6i\y,mi) (xp(y\6 i ,m i )p{e i \mi), (2) 
and thus results from integrating the un-normalised posterior across the 0j parameter space, 

p(y\rrii)= p(y\0i,mi)p(0i\mi) dOi. (3) 



This of course assumes that the prior distribution for 6i is proper. The marginal likelihood is often then used 
to calculate Bayes factors when one wants to compare two competing models, m; and nij, 



_ p(y\m) _ p(mj\y) pijnj) 
p(y\mj) p(mj\y) p(m) 



BFa = T 1 = *) 3 ' (4) 



Here, p(rrii\y) is the posterior probability for model m; and it can be evaluated, using the evidence for each 
of the collection of models under consideration, 

p(mi\y) oc p(y\mi)p(mi), for i = 1, . . . , I (5) 

Estimation of the evidence is a non-trivial task for most statistical models and there has been consider- 
able effort in the literature to find algorithms and methods for this purpose. Laplace's method (Tierney and 
Kadane, 1986) is an early approach and very widely used. Other notable and popular approaches include 
Chib's method (Chib 1995), annealed importance sampling (Neal 2001), nested sampling (Skilling 2006), 
bridge sampling (Meng and Wong, 1996) and power posteriors (Friel and Pettitt, 1998) which is the focus 
of this paper. For a recent review and perspective on these and other methods, see Friel and Wyse (2012). 

This paper is organised as follows. Section |2] outlines the power posterior method, and the approach we 
propose to improve estimation of the evidence. Section [3] illustrates the potential gain from implementing 
the methodology which we propose. We offer some conclusions in Section [4] 



2 The power posterior approach 

In what follows we will drop the explicit conditioning on model rm for notational simplicity. We follow the 
notation of Friel and Pettitt (2008) and denote the power posterior by 

p t (6\y) oc p(y\efp(e), t e [0,1] (6) 

with%|i) = f pivlOfpffldB. (7) 
Je 

where t 6 [0, 1] is thought of as a temperature, which has the effect of tempering the likelihood, whereby 
at the extreme ends of the temperature range, po(9\y) and pi(6\y) correspond to the prior and posterior, 

2 



respectively. The power posterior estimator for the evidence relies on noting that 

r ± P (y\ey P (6)d6 



z(y\t) Je dt 
1 



p(y\6) t log(p(y\e))p(9)d9 

zm log { p(y\9))d9 
E \ y , t ]og(p(y\O)). (8) 



z(y\t) Je 

p{y\6fp{e) 



As a result 

Be\ y , t log(jp{y\e))dt = [\og(z{y\t))}l 

= \og{z{y\t = 1)) (assuming that the prior is normalised) (9) 

which is the log of the desired marginal likelihood. 

In practice the temperature range is discretised as = to < h, ■ ■ ■ , t n = 1 to form an estimator based on 
(©. For each t,- L , a sample from p(6\y, U) can be used to estimate Egu^. log(p(y\9)). Finally, a trapezoidal 
rule is used to approximate 

iog P « - f> - *,-0 ^ Mm \ + E<1 "" '° gW#)> ) ■ 00) 

Discretising i introduces an approximation into this method and the two goals of this paper are to reduce the 
bias in the power posterior estimation method due to the approximation and also to find an adaptive method 
for choosing the temperature rungs required. For both of these we will exploit the fact that the gradient of 
the expected log deviance curve equals its variance, as we now outline. 
Differentiating ^e\y,t log(p(y|0)) with respect to t yields 

d -E fl |„ it log(p(y|0)) = j\og{p{y\8))^p t {9\y)d8 

1 Pt(0\y)d8 
Pt(0\y)d6 



dt °™ " Je °™ "dt 

log(p(#)) [log(p(y|0)) - -±^±z(y\t) 



log(p(y|0)) \log(p(y\6)) - ^ \og(z(y\t)) 
I dt 

= E %ji log(p(y\9)) 2 - (E % , t \og(p(y \9))) 2 

= V elyjt (log(p(y\6))) (11) 
where V e i y t (log(p(y\9))) denotes the variance of the log deviance at temperature t. 
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2.1 Reducing the bias by improving the numerical integration 

Equation (fTTT t immediately provides two useful pieces of information. First, the curve which we wish 
to integrate numerically is (strictly) increasing. Secondly, we can improve upon the standard trapezium 
rule used to numerically integrate the expected log deviance by incorporating derivative information at 
virtually no extra computational cost (the cost merely of calculating the variance of a set of simulations 
for fixed t). We do this by using the corrected trapezium rule which comes from an error analysis of the 
standard trapezium rule, see for example Atkinson and Han (2004), Section 5.2; when integrating a function 
/ between points a and b 



f(x)dx = (b — a) 



f(b) + f(a) 



12 



-f"(c) 



(12) 



where c is some point in [a,b]. The first term of the right hand side of this equation is the usual trapezium 
rule and the second can be approximated using 

/'(*) " /'(«) 



so that 



f"(c) 
f(x)dx 



(b-a) 



f(P) + f(a) 
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[/'(*) " /'(«)] • 



(13) 



This latter form motivates the corrected trapezium rule which for unequally spaced x-axis points, taken 
together with the information derived above regarding the derivative of the log deviance gives 

"E %A log(p(y\9)) + E %A+1 log(p(y|0))" 



log(z(y\t = 1)) 



n-1 

(u+i - u) 

i=0 

1 - 1-* - 



E 

i=0 
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V % , ti+1 (log(p(2/|0))) - V %)tt (log(p(l/|0))) (14) 



where both the expectations {Eg\ yit . log(p(y\9))} and variances {V d \ yt .(log(p(y\9)))} are to be estimated 
using MCMC runs at a number of values of ij. 



2.2 Adaptive choice of the temperature placement 

The next question which arises is how to choose the {ij} between to = and t n = 1. Friel and Petti tt (2008) 
find that setting ti = (i/n) 5 performs well. We refer to this as the powered fraction (PF) schedule. Lartillot 
and Philippe (2006) discuss very similar ideas in the phylogenetics literature, although using Simpson's rule 
for the numerical integration; they use equally spaced temperatures between and 1 . 

Here we will only consider the discretisation error associated with the numerical integration, rather 
than the stochastic error arising with sampling from the different p t . (6\y). Calderhead and Girolami (2009) 
show that this discretisation error depends upon the Kullback-Liebler distance between successive p ti (9\y). 
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Lefebvre, Steele and Vandal (2010) also consider a symmetrised Kullback-Liebler divergence in picking op- 
timal schedules for path sampling. At first glance the Kullback-Liebler distance does not seem a particularly 
tractable quantity to manipulate. However, these papers and Behrens, Friel and Hum (2012) all note that, in 
the notation of this paper, 

n-l 

J2(KL\p t My),Pt i+1 (o\y)] + KL\p ti+1 (e\y),PtMy)}) = 2S„(t , ■■■,*«) (15) 

i=l 

where KL denotes the Kullback-Liebler distance and 

n— 1 n— 1 

S n (t , ...,t n ) = 2(<<+i " **)E%,t j+1 log(p(y|0)) - ~ ^) E %,t« lo g(K#))- ( 16 ) 

i=0 i=0 

iS n can be interpreted graphically as the sum of the rectangular areas between a lower and an upper approxi- 
mation to the integral of ^e\y,u l°g(p(y 1$)) between to = and ii = 1. Behrens, Friel and Hum (2012) use 
minimising S n as a rationale for choosing the temperatures in tempered transitions. We propose to use the 
same target in selecting the {ti} for power posteriors. However, unlike in tempered transitions where the 
tuning forms a small part of the overall computational load, here the cost is almost exclusively the estima- 
tion of ^e\y,u l°g(p(y|$))- We propose the following scheme: Initialise a set of m {ti} using the geometric 
placing including and 1 (we will see in later examples why a reasonable starting point is necessary) where 
m is a small proportion of the proposed total number of rungs n. These m {ti} contribute (m — 1) terms 
{[t i+ i - ti]\E e \ ytt . +1 \og(p(y\0)) - Eo\ y ,ti log(p(y|0))]} which sum to give S n . Identify the largest of these 
terms and locate the next point in the corresponding interval, say [t k ,t k+ i]. Since we do not want to use 
computational resources in performing a search for the optimal location of the new ti (there is no analytic 
solution), we follow a low cost route using the estimated gradients/variances at t k and t k+ i. If the estimated 
gradient at t k is denoted by V*. and that at t k+ \ by Vk+i, we set the new point to be 

t = t k + Vk +} r (t k+1 - t k ). (17) 
Vfc + V k+l 

This scheme will almost certainly not identify the optimal placing of the n rungs. However it is quick, cheap 
and intuitively reasonable. (In practice, Monte Carlo error can mean that the function is not increasing and 
so the criterion is changed to picking the interval with the largest absolute contribution to S n .) 



3 Examples 

We present three examples which illustrate the gains that arise from employing the methods developed here. 
The first example is a non-nested linear regression comparison for which the marginal likelihoods can be 
calculated analytically. Example 2 is a larger problem, choosing between two logistic regression models, for 
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which an analytic solution is not possible. These first two examples were included in the review paper by 
Friel and Wyse (2012) where the performance of power posteriors was compared to other existing methods. 
The final example is by far the largest and exhibits the most interestingly shaped ~Eg\ yj t log(p(y|0)). 

3.1 Example 1: Radiata pine 

The first example compares two linear regression models for the Radiata pine data originally in Williams 
(1959). The response variable here is the maximum compression strength parallel to the grain, j/j, while the 
predictors are density, %{, or density adjusted for resin content, z%, for n = 42 specimens of radiata pine. 
Two possible Gaussian linear regression models are considered; 

Model 1: y; t = a + (3(xi - x) + Cj, e; ~ iV(0, r _1 ), i = l,...,n, 
Model 2: m = 7 + 6(2% - z) + rji, rji ~ N(0, X^ 1 ), i = l,...,n. 

Priors are chosen to match the analyses of Friel and Wyse (2012) (baring a notational factor of 2). The 
regression parameters (a, (3) T and (7, 5) T are taken to be Normally distributed with mean (3000, 185) T and 
precision tQq and XQq respectively where Qq = diag(ro, so)- The values of ro and so were fixed to be 0.06 
and 6. A gamma prior with shape ao = 3 and rate bo = 2 x 300 2 was assumed for both r and A. 

Following the comparisons of Friel and Wyse (2012), we consider estimating the evidence using 10, 
20, 50, 100 or 200 rungs in the tempering scheme. The parameters at all levels are updated using the 
Gibbs sampler. For this example, both the adaptive and the PF spacings use 20000 iterations at each rung, 
discarding the first fifth of these as burn in. Figure [TJ shows the expected log deviance curves for the two 
models using 200 rungs, their shapes suggesting that PF spacing might perform competitively (Behrens, 
Friel and Hum (2012) show that a scheme where ti/ti+i is a constant for i > minimises S n when the 
integrand takes the form ^ + K2 for some constants K\ and K2). 

Figure |2] shows the upper and lower bounds of the evidence (in black), the uncollected estimate (in 
red) and the corrected estimate (in blue) all for model 1 as the number of rungs increases. The PF spacing 
results are denoted by solid lines and the adaptive spacing results by dashed lines. The true value of the 
evidence is known for this example and is marked by a horizontal line. As the vertical scale differs quite 
significantly between n = 10 and n = 200, the figure is split into two plots, small numbers of rungs (where 
the upper and lower bounds are not tight) and large numbers of rungs. The adaptive temperature placement 
is initialised using the 10 rung PF placement. Since for the adaptive spacing, increasing the number of rungs 
by one requires only one additional set of MCMC iterations at the new temperature, there is an averaging 
effect and the dashed lines appear smoother than the solid ones (where for an increase of one rung, all 
temperatures apart from to = and t n = 1 change and so the estimates at successive rungs are independent 



6 



The estimated expected log deviance for both models 
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Figure 1 : The expected log deviance curves for the two Radiata models using 200 rungs. 
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Figure 2: For the Radiata example. Upper and lower bounds (in black), uncorrected estimates (red), cor- 
rected estimates (blue) as the number of rungs increases for model 1. Solid lines indicate PF spacing, dashed 
lines the adaptive schedule. 8 
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Table 1: Estimated bias (and standard deviation) in estimating the evidence for the two Radiata pine models. 

of one another). From this figure, it appears that the corrected estimates converge faster towards the true 
value than do the uncorrected ones initially. By construction, the adaptive and PF schedules coincide at 10 
rungs. Immediately after that, the adaptive schedule initially provides wider bounds on the evidence but 
after approximately 25 rungs, the bounds are consistently narrower. 

To quantify these observations, the bias is estimated as per the approach of Friel and Wyse (2012), 
performing 50 replicates at 10, 20, 50, 100 and 200 rungs using 10000 iterations of which the first fifth are 
discarded as burn in. The average and standard deviation of the 50 biases are given in Table [T] As might be 
expected from the concave shape of the log deviance curves, the uncorrected integrals tend to underestimate 
the evidence, giving negative biases which decrease as the number of rungs increases. To visualise the 
effect of both the correction and the adaptive placing of temperatures, Figure [3]plots the 50 observed biases 
separately for each number of rungs and under the two spacings with or without correction. Correction is 
particularly effective when smaller numbers of rungs are being used. The final panel illustrates the effect of 
the correction when using 100 rungs by plotting the corrected points against their uncorrected values (the 
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line y = x is shown in red, corresponding to no effect of correction). There is an interesting difference 
between the adaptive and the PF versions, less correction is needed for the adaptive schedule, that is, it is 
doing a better job of the numerical integration. 

Given the good reductions in bias seen in Table [TJ it is important to ask how much extra time is required. 
To assess this, a total of 10 runs for model 1 using 20000 iterations and 200 temperatures were timed. Four 
versions of the algorithm were considered, corresponding to Table [T] All the coding was in R and times are 
given relative to the PF non-corrected version: 

PF non-corrected PF corrected Adaptive non-corrected Adaptive corrected 
1.0000 1.0054 1.0017 0.9986 

The adaptive selection of temperatures and the correction term in the numerical integration come at negli- 
gible computational cost. Given the reductions in bias achievable by the correction in particular, there is no 
reason at all not to adopt this modification. 

3.2 Example 2: Pima indians 

We turn next to the Pima Indian example considered by Friel and Wyse (2012), originally described by 
Smith et al (1988) . These data record diabetes incidence and possible disease indicators for n = 532 Pima 
Indian women aged over 20. The seven possible disease indicators are the number of pregnancies (NP), 
plasma glucose concentration (PGC), diastolic blood pressure (BP), triceps skin fold thickness (TST), body 
mass index (BMI), diabetes pedigree function (DP) and age (AGE), with all these covariates standardised. 
The model assumed for the observed diabetes incidence, y = (yi, . . . , y n ), is 

n 

p(y\o) = l[pT( 1 -Pi) 1 ~ y ' (18) 

i=i 

where pi is the probability of incidence for person i, and p. t is related to the i th person's covariates and a 
constant term, denoted by Xj = (1, xn, . . . , Xid) T , and the parameters, 9 = (9 , 9\, . . . , 9d) T , by 

log (t3^) = ° Tx i (I 9 ) 

where d is the number of explanatory variables. An independent multivariate Gaussian prior is assumed for 
9, with mean zero and non-informative precision of r = 0.01, so that 

p(9) = (27r)- d /V/ 2 exp • (20) 

There are 129 potential models (2 7 models with covariates plus a model with only a constant term). A 
long reversible jump run (Green, 1995) revealed the two models with the highest posterior probability: 
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Figure 3: The 50 observed biases for model 1 using different numbers of rungs and the four schemes for 
the Radiata data: l=Uncorrected PF spacing; 2=Corrected PF spacing; 3=Uncorrected adaptive spacing; 
4=Corrected adaptive spacing. The red line in the comparison plot corresponds to the correction having no 
effect. 



Model 1: logit(p)= 1+NP+PGC+BMI+DP 
Model 2: logit(p)= 1+NP+PGC+BMI+DP+AGE 

Figure [4] shows the estimated log deviance curves for these two competing models. Although the shapes 
are roughly similar to those in the Radiata example, notice the difference in scale on the y-axis compared to 
FigureQ] For these two models, the power posterior is not amenable to the Gibbs sampler and so a Metropolis 
update is used instead. This raises the problem of proposal scaling at the different temperatures. Since both 
the correction and the adaptive temperature placements rely on having good estimates of the variance of the 
log deviance, mixing is an important issue here. As an alternative to the approach taken in Friel and Wyse 
(2012), we have considered here a joint update of all the model parameters using a multivariate Normal 
proposal centred at the current value and with diagonal variance matrix, entries min(0.01/i, 1/r) where t 
is the temperature and r is the precision of the prior. Run lengths are doubled compared to the previous 
example, using 40000 iterations and discarding the first fifth as burn in, with the adaptive temperature 
placement still initialised at the 10 rung PF placement. 

Figure[5]shows the comparison of a single adaptive run with the corresponding PF scheme. Unlike in the 
Radiata example where the evidence can be evaluated analytically, Friel and Wyse (2012) use the Laplace 
approximation of the log evidence (Tierney and Kadane, 1986) as the "benchmark" in assessing convergence 
and bias. However this is not necessarily very accurate and so we replace the Laplace approximation by 
a very long run (2000 rungs) of the power posterior approach; the estimates we get are —257.2342 and 
—259.8519 for models 1 and 2 respectively as opposed to the Laplace approximations of —257.2588 and 
—259.8906. Comparing Figures |2] and \5\ the greater y-axis scale of the expected log deviance curve for 
this example has led, not surprisingly, to wider intervals for the evidence. As before though, the adaptive 
temperature scheme initially under-performs the PF scheme but as the number of rungs increases beyond 
about 25 it generates narrower upper and lower bounds. The corrected estimates also appeal - to converge 
faster than the uncorrected ones towards the benchmark log evidence. 

Table [2] shows the estimated biases and standard deviations (here estimated using 25 replicates rather 
than the 50 of the previous example for reasons of speed). These estimated biases are also depicted in 
Figure [6] Both the correction and the adaptation seem effective in reducing bias, the former particularly 
dramatically for small numbers of rungs while for the latter there is the same "catch-up" effect comparing 
the uncorrected versions for small numbers of rungs but thereafter there is a clear benefit. Interestingly, as in 
the previous example, the correction evens out the differences between the PF and the adaptive schedules' 
biases, although as the last panel in Figure [6] shows, to do so requires a larger correction for the PF. 
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The estimated expected log deviance for both models 
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Figure 4: The expected log deviance curves for the two Pima models using 200 rungs. 
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Results using 10 to 100 rungs 





Figure 5: For the Pima example. Upper and lower bounds (in black), uncorrected estimates (red), corrected 
estimates (blue) as the number of rungs increases for model 1 . Solid lines indicate PF spacing, dashed lines 
the adaptive schedule. 14 
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Figure 6: The 25 observed biases for model 1 using different numbers of rungs and the four schemes for the 
Pima Indian data: l=Uncorrected PF spacing; 2=Corrected PF spacing; 3=Uncorrected adaptive spacing; 
4=Corrected adaptive spacing.The red line in the cojnparison plot corresponds to the correction having no 
effect. 
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(0.1640) 


/ r\ Anoi\ 

(0.0782) 


(0.0706) 


(0.0386) 


(0.0331) 




Adaptive corrected 


U.OOJO 


U.Uo3o 


a m^o 


A A1 OA 


A AAQO 






(0.1398) 


(0.0791) 


(0.0703) 


/ r\ r\n O A \ 

(0.0384) 


(0.0331) 


Model 2 


PF uncorrected 


A a Ann 

-4.14/ / 


-1.U254 


A 1 £ /I O 

-U. 1643 


A A/l 

-U.U433 


A A1 A A 

U.U14U 






(0.2028) 


(0.1172) 


(0.0726) 


(0.0359) 


(0.0466) 




PF corrected 


A A^'2'2 


A AC] O 

U.U538 


A AAA A 
0.0004 


-U.UU25 


A AA'SA 

-0.003y 






(0.1697) 


(0.1074) 


(0.0715) 


(0.0357) 


(0.0465) 




Adaptive uncorrected 


-4.1619 


-1.3083 


-0.1165 


-0.0281 


-0.0020 






(0.1606) 


(0.1025) 


(0.0521) 


(0.0458) 


(0.0316) 




Adaptive corrected 


1.0125 


0.1045 


-0.0057 


-0.0020 


0.0041 






(0.1504) 


(0.0980) 


(0.0515) 


(0.0457) 


(0.0316) 



Table 2: Estimated bias (and standard deviation) in estimating the evidence for the two Pima Indian models. 



3.3 Example 3: Galaxy data 

To demonstrate a large application with a more challenging integral than the previous two, we use the much- 
studied Galaxy data set, see for example Richardson and Green (1997), which comprises measurements on 
the velocities of 82 galaxies. Denoting the 82 measurements by y = {y\, . . . , ys2}, we follow Richardson 
and Green (1997) in incorporating corresponding latent allocation variables z = {z\, . . . , z^}- Given 
z i = J. Hi follows the j th of the k component Gaussian distributions of the mixture, 

p(yi\zi = j, fij, a]) = -jL= exp ( ) * = 1, . . . , 82. (21) 

Conditional independence is assumed for the {yi} and we specify independent standard proper priors: 

k 

P( z i = j) = w ji where wj = 1 (22) 

i=i 

{wi,...,Wk} ~ Dirichlet(l,k) (23) 
N ~ iV(0,1000), j = l,...,k (24) 
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o a ~ InvGam(l,l), j = l,...,k. (25) 

The weights, means and variances are all updated using the Gibbs sampler but we use a Metropolis algorithm 
with a discrete uniform proposal for the allocation variables. Run lengths of 40000 will be used, discarding 
the first tenth as burn in. Behrens, Friel and Hum (2012) considered this model when studying temperature 
placement for the tempered transition algorithm, finding that its expected log deviance curve had some 
interesting features. Figure |7]shows the estimated log deviance curves for k = 1 through to k = 7. The upper 
panel shows just the k = 1 curve as they are virtually indistinguishable due to the huge scale change as t 
approaches zero. However by restricting the y-axis, interesting differences can be seen for larger temperature 
values in the lower panel. Capturing these features efficiently for the numerical integration as k changes is 
the challenge. The top panel of Figure [8] shows how the adaptive scheme differs from the PF one for 200 
rungs using k = 1 and k = 3. For k = 1, the adaptive scheme moves the points more densely towards zero 
where the expected log deviance decreases very rapidly. For k = 3, the expected log deviance also has a 
section of rapid change in the mid-range temperatures and the adaptive scheme attempts to capture this. 

The massive gradient of the expected log deviance curves illustrates a potential pitfall of the adaptive 
temperature scheme. Figure [9] again shows a comparison between a rung of adaptive placement against 
the uncorrected PF scheme as the numbers of rungs increases. In the upper figure, the adaptive scheme 
is initialised at 10 PF spacings and in the lower at 20 PF spacings. In both cases by the time 200 rungs 
are being used, the adaptive version outperforms the PF in terms of the separation of the upper and lower 
bounds. But there is a dramatic difference at the start; the 10 rung initialisation takes far longer to catch 
up. At the end of the process, the two sets of {U} are not very different (Figured bottom panel). However 
the rate of change of the curve towards t = is so great that the 10 rung version starts by filling in lots of 
temperatures close to zero, with the large differences in estimated gradients there meaning that our location 
procedure for additional points locates them very close together initially. There is at least an immediate 
diagnostic that the initial number of rungs is insufficient in that the corrected estimate lies outside the lower 
and upper bounds (Figure |9j upper panel); the changes in the gradient are so huge between t = and the 
next few {ti} that pairwise differences do not form a good estimate of the second derivative (Equation (fT3T>). 
This effect is the reason why we do not initialise the process using just the t = and t = 1 points which are 
common to all schemes. 

Using a 20 rung initialisation with 200 rungs in total, Figure \10\ shows results for the mixture models 
with k = 1 to k = 7. The highest log evidence belongs to the k = 3 model, shortly followed by k = 4, 5, 6 
and 7 in that order. The first set of non-overlapping discretisation bounds is between k = 3 and k = 7; the 
difference here between the corrected estimates is 3.38, with a corresponding Bayes factor of 29.48 for this 
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The estimated expected log deviance 
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The estimated expected log deviance, restricted y-axis 
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Figure 7: The expected log deviance curves for the Galaxy models using 200 rungs. 
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Temperature placements for the Galaxy mixture 
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Figure 8: Top panel: The placement of the 200 rung schedules for the Galaxy data; the straight line indicates 
an exact match between PF and adaptive schemes. Bottom panel: The placements for k = 3 changing the 
initialisation from 10 to 20 rungs. 19 



Results using 10 rung initialisation 
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Results using 20 rung initialisation 
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Figure 9: Upper and lower bounds (in black), uncorrected estimates (red), corrected estimates (blue) as the 
number of rungs increases for the k = 3 mixture model. Solid lines indicate PF spacing, dashed lines the 
adaptive schedule. 20 



choice of priors. The adaptive scheme has the additional benefit that if these discretisation bounds are still 
considered too wide after using the anticipated maximum number of temperatures n, the process can simply 
be run on with additional temperatures placed by the same algorithm. The same cannot be said of the PF 
scheme where it is not immediately clear how to add additional temperatures. 

4 Conclusions 

This article has, we hope, illustrated the potential gains that can be made when estimating the evidence using 
power posteriors by correcting the numerical integration error and by adaptively choosing the temperature 
ladder. The methods that we have outlined come at virtually no extra computational cost, and we would 
therefore recommend that these are routinely used when implementing the power posterior approach. Given 
how effective the correction is but remembering that it does rely on good estimates of variance, it may be 
better to use a moderate size of n with long MCMC runs at each U, rather than dividing up the same total 
number of iterations into short runs with a large n. 

What this article does not do is to give guidance as to how to allocate computational resources between 
the different temperatures. We have seen in our examples that the gradient and thus also the variance of 
~E>o\ y t t log(p(y|#)) i s largest as t — > 0, suggesting that we should allocate more MCMC iterations here rather 
than as t —> to get good estimates. On the other hand, when t is small, the power posteriors pt{0\y) 
will probably be easy to sample compared to when t = 1 and so we should also take into account the 
MCMC effective samples sizes. Neither the gradients nor the effective sample sizes can be known before 
sampling is earned out! There probably is some scope for an adaptation scheme here too, perhaps allocating 
some fraction of the total number of MCMC iterations evenly over the temperatures before allocating the 
remainder based on what we have learned in this initial phase. This point also reinforces our caveat: what 
we have addressed here is discretisation error, the bounds we give are (noisy) bounds on this error and not 
credible intervals in the usual sense. 

In our examples, applying the correction term has effectively smoothed over the benefits of the adaptive 
scheme over the PF one (just working a little harder in the latter case). This numerical analysis trick is 
peculiar to this particular use of tempered distributions. In general we suspect though that the adaptation 
ideas developed here could find wider use in other tempered schemes described in the literature where 
numerical integration is not involved. 
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The estimated log evidence for varying k 
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Figure 10: The estimated log evidence for the Galaxy data. Top panel: Dashed lines indicate upper and 
lower bounds, solid lines indicated the corrected estimate of the log evidence for k = 1 through to k = 4. 
Bottom panel: vertical bars indicate the width of the i®ver and upper bound, circles the corrected estimates, 
k = 2 through to k = 7. 
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